//------------------------------------------------------------------------------
//!	@file	dsp10.c
//!	@brief
//!		10-Bit DSP processing functions
//!
//!	@Description
//!		10-Bit DSP processing functions
//!
//!
//!	Originally Developed: July, 1997 to October, 1999
//!
//!	$Id: dsp10.c 158 2013-03-03 07:00:15Z ac.verbeck@gmail.com $
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//	Copyright and Disclaimer Notice
//
//	Copyright (c) 2013 A.C. Verbeck
//
//	This file is part of CreateLUT: A LUT creation library.
//
//	CreateLUT is free software: you can redistribute it and/or modify
//	it under the terms of the GNU General Public License as published by
//	the Free Software Foundation, either version 3 of the License, or
//	(at your option) any later version.
//
//	CreateLUT is distributed in the hope that it will be useful,
//	but WITHOUT ANY WARRANTY; without even the implied warranty of
//	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//	GNU General Public License for more details.
//
//	You should have received a copy of the GNU General Public License
//	along with CreateLUT.  If not, see <http://www.gnu.org/licenses/>.
//
//------------------------------------------------------------------------------

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <memory.h>

#include <stdint.h>
#include <stdbool.h>

#include "dda10.h"
#include "SysParam.h"

#include "dsp10.h"

typedef struct FLAT_SPOTtag {
	int32_t Start;
	int32_t	Rise;
	int32_t	Run;
} FLAT_SPOT;

extern bool NoDSP;

//
//	Name:
//		DSP10_BuildBitLUT
//
//	Description:
//		Generates a single LUT based on the following:
//		-- A normaized LUT
//		-- Gain
//		-- Gamma
//
//	Parameters:
//		uint8_t*	OutLut:		Output LUT (256 entries)
//		double*	NormLut:	Input Normalized LUT (256 entries)
//		double	Gain:		Gain to be applied to this LUT
//		double	Gamma:		Gamma to be applied to this LUT
//
//	Returns:
//		(No direct return)
//		Indirect: The OutLut[TableSz] array is filled.
//
void DSP10_BuildLUT(uint8_t* OutLut, double* NormLut, double Gain, double Gamma)
{
	uint16_t	t1[CURVE_SZ_10BIT];
	double	tIdx,tInc;
	double	dVal;
	int32_t	i;
	int32_t	sml,lrg;

//
//	Pre LUT production
//
	tInc = (double)1/(double)(CURVE_SZ_10BIT-1);

//
//	LUT production (produce a 10bit LUT)
//
	for (i=0,tIdx=0.0; i<CURVE_SZ_10BIT; i++,tIdx+=tInc) {
		dVal = pow(tIdx,Gamma)*Gain;
		t1[i] = DSP10_NormalLookUp(NormLut, dVal);
	}

//
//	Send to output LUT
//
	for (sml=0,lrg=0; sml<CURVE_SZ_08BIT; sml++,lrg+=4) {
		if (sml<128)
			OutLut[sml] = (uint8_t)t1[lrg];
		else
			OutLut[sml] = (uint8_t)(t1[lrg]/4);
	}
}


//
//	Name:
//		DSP10_Normalize
//
//	Description:
//		Normalize the InLut to OutLut.
//		Locates the Min and Max for normalization.
//		Applies Min and Max to each value to range
//		it between 0..1.
//
//	Parameters:
//		double*	OutLut:		Output Normalized LUT
//		double*	InLut:		Input LUT
//
//	Returns:
//		(No direct return)
//		Indirect: The OutLut[TABLE_SZ] array is filled.
//
void DSP10_Normalize(double* OutLut, double* InLut)
{
	double	MinVal,MaxVal;
	double	delta;

	int32_t	i,j;

	MinVal = InLut[0];
	MaxVal = InLut[0];

//
//	Locate Min & Max Value
//	Note: input table is 8 bits
//
	for (i=0; i<CURVE_SZ_08BIT; i++) {
		if (MinVal > InLut[i]) MinVal = InLut[i];
		if (MaxVal < InLut[i]) MaxVal = InLut[i];
	}

//
//	Normalize as an 8 bit extended array (fill every 4th value)
//	Input array is a 256 element array
//
	for (i=0,j=0; i<CURVE_SZ_08BIT; i++,j+=4)
		OutLut[j] = (InLut[i]-MinVal) / (MaxVal-MinVal);

//
//	Fill in the missing values using linear interpolation
//	Do all but last one.
//
	for (i=0; i<CURVE_SZ_10BIT-4; i+=4) {
		delta = (OutLut[i+4]-OutLut[i])/4;
		for (j=i+1; j<i+4; j++) {
			OutLut[j] = OutLut[j-1]+delta;
		}
	}

//
//	Interpolate last three points
//
	delta = (1.0-OutLut[i])/4;
	for (j=i+1; j<i+4; j++) {
		OutLut[j] = OutLut[j-1]+delta;
	}
}


//
//	Name:
//		DSP10_Filter00
//
//	Description:
//		Filter the function by sorting
//
//	Parameters:
//		double*	OutLut:	Destination array for data (256 entries)
//		double*	InLut:	Source array for data (256 entries)
//
//	Returns:
//		(none)
//
int dbl_cmp(const void* a, const void* b)
{
	return ((*(double*)a > *(double*)b) ? 1:-1);
}
void DSP10_Filter00(double* OutLut, double* InLut)
{
	memcpy(OutLut,InLut,CURVE_SZ_10BIT*sizeof(double));
	if (NoDSP == true) return;
	qsort((void*)OutLut, CURVE_SZ_10BIT, sizeof(double), dbl_cmp);
}


//
//	Name:
//		DSP10_NormalLookUp
//
//	Description:
//		Locates the closest LUT index in the normalized LUT.
//		This is done by locating the closest LUT value to input value
//		and returning the index.
//
//	Parameters:
//		double*	NormLut:	Normalized LUT
//		double	InVal:		Input value to be located
//
//	Returns:
//		Index of the closest value.
//
uint16_t DSP10_NormalLookUp(double* NormLut, double InVal)
{
	double	MinDelta;
	int32_t	MinIdx,i;

	MinIdx = 0;
	MinDelta = fabs(NormLut[0]-InVal);

	for (i=0; i<CURVE_SZ_10BIT; i++) {
		if (MinDelta > fabs(NormLut[i]-InVal)) {
			MinIdx = i;
			MinDelta = fabs(NormLut[i]-InVal);
		}
	}
	return (uint16_t)MinIdx;
}


//
//	Name:
//		DSP10_Smooth
//
//	Description:
//		Uses a single octant dda to do a linear smooth of any
//		flat spots on InLut.
//
//	Parameters:
//		OutLut[CURVE_SZ_10BIT]:	Output LUT
//		InLut[CURVE_SZ_10BIT]:	Input LUT
//
//	Returns:
//		(No direct return)
//		Indirect: The OutLut[CURVE_SZ_10BIT] array is filled.
//
void DSP10_Smooth(uint16_t* OutLut, uint16_t* InLut, int32_t Threshold)
{
	int32_t	i,fIdx;
	int32_t	Dy[CURVE_SZ_10BIT];
	FLAT_SPOT Fs[CURVE_SZ_10BIT];

	memcpy(OutLut,InLut,CURVE_SZ_10BIT*sizeof(uint16_t));
	if (NoDSP == true) return;

//
//	Initialize
//
	for (i=0; i<CURVE_SZ_10BIT; i++) {
		Fs[i].Start = 0;
		Fs[i].Rise = 0;
		Fs[i].Run = 0;
	}

//
//	Create the difference array.
//
	for (i=0; i<(CURVE_SZ_10BIT-1); i++) Dy[i] = InLut[i+1]-InLut[i];

//
//	Find the flat spots
//
	fIdx=0;
	for (i=0; i<CURVE_SZ_10BIT-1; i++) {
		if (Dy[i] == 0) {							// If we are at a flat spot
			if (Fs[fIdx].Run == 0) {				// And this is the first time we have hit this spot
				Fs[fIdx].Start = i;					// Set start to first point32_t of flat spot
				Fs[fIdx].Rise = InLut[i];			// Save the current point32_t
				Fs[fIdx].Run = 1;					// Set length to 1
			} else {									// More than one point32_t on this flat spot
				Fs[fIdx].Run++;
			}
		} else {													// Not at a flat spot
			if (Fs[fIdx].Run != 0) {
				Fs[fIdx].Rise = (int32_t)InLut[i+1]-Fs[fIdx].Rise;	// Get the rise
				fIdx++;											// Move to next entry (if we need to)
			}
		}
	}

//
//	Smooth output
//
	for (i=0; i<fIdx; i++) {
		if (Fs[i].Run > Threshold) {
			int32_t x1 = Fs[i].Start;
			int32_t y1 = OutLut[x1];
			int32_t x2 = x1+Fs[i].Run;
			int32_t y2 = y1+Fs[i].Rise-1;
			dda10(OutLut, x1, y1, x2, y2);
		}
	}
}


//
//	Name:
//		DSP10_Curve
//
//	Description:
//		Does a parabola curve fit on the bottom part of each LUT.
//
//	Parameters:
//		OutLut[CURVE_SZ_10BIT]:	Output LUT
//		InLut[CURVE_SZ_10BIT]:		Input LUT
//		Threshold:				The end point of the curve fitting.
//		Shape:					This is the shape of the parabola
//								(Good values range from 1.5 to 2.0)
//
//	Returns:
//		(No direct return)
//		Indirect: The OutLut[CURVE_SZ_10BIT] array is filled.
//
void DSP10_Curve(uint16_t* OutLut, uint16_t* InLut, int32_t Threshold, double Shape)
{
	int32_t	i, X_Th;
	double	t1;

	memcpy(OutLut,InLut,CURVE_SZ_10BIT*sizeof(uint16_t));
	if (NoDSP == true) return;

//
//	Locate the threshold in the LUT
//
	if (Threshold > InLut[CURVE_SZ_10BIT-1]) {
		printf("DSP10_Curve: Warning: threshold too high: %d\n",Threshold);
		return;
	}

	for (i=0; i<CURVE_SZ_10BIT; i++) {
		if (InLut[i] > Threshold) {
			X_Th = i;
			break;
		}
	}

//
//	Do the curve fit
//
	for (i=0; i<X_Th; i++) {
		t1 = pow(((double)i/(double)X_Th), Shape) * (double)Threshold;
		t1 += 0.5;
		OutLut[i] = (uint8_t)t1;
	}
}

//
//	End: dsp10.c
//